# Temas para os gráficos
theme.base <- theme_minimal(base_size = 11) +
theme(
axis.text = element_text(size = 8),
plot.title = element_text(hjust = 0.5, size = 11, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, size = 9),
plot.caption = element_text(size = 8),
axis.title = element_text(size = 8),
legend.title = element_text(size = 8),
panel.grid.major = element_line(colour = "grey90", linewidth = 0.5),
panel.grid.minor = element_line(colour = adjustcolor("grey90", alpha.f = 0.4), linewidth = 0.5),
panel.border = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
axis.line.x = element_line(colour = "grey"),
axis.line.y = element_line(colour = "grey"),
)
theme.no_legend <- theme(legend.position = "none")
theme.no_grid <- theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
theme.no_axis <- theme(
axis.line.x = element_blank(),
axis.line.y = element_blank()
)
theme.no_axis_title <- theme(
axis.title.x = element_blank(),
axis.title.y = element_blank()
)
# Theme for timeseries
theme.ts <- theme.base +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()
) +
theme.no_legend
plotly.base <- function(p) {
p %>%
layout(
margin = list(b = 60, t = 80)
) %>%
config(mathjax = 'cdn')
}# Função para salvar e carregar resultados
cache_dados <- function(chave, funcao_geradora) {
c <- paste(chave, "rds", sep = ".")
cache <- file.exists(c)
if (!cache) {
resultado <- funcao_geradora()
saveRDS(resultado, c)
} else {
resultado <- readRDS(c)
}
resultado
}βAR(1) é um modelo de séries temporais autorregressiva de ordem 1 com distribuição Beta – onde a variável de interesse está restrita no intervalo (0, 1) –, este modelo foi proposto por Rocha e Cribari-Neto:
O objetivo do presente trabalho é aplicar métodos de controle de processos βAR(1) para detectar mudanças no processo. Para isso, utilizamos o fato de que, quando modelada corretamente, os resíduos de uma série temporal, como a βAR(1), são independentes e normalmente distribuídos com média zero e variância constante.
Desta forma, quando o processo sofre uma mudança, espera-se que os resíduos não sejam mais independentes e que a média e a variância dos resíduos sejam diferentes. Assim, podemos utilizar métodos de controle de processos para detectar essas mudanças.
Olhando sob a perspectiva de teste de hipóteses, temos:
Além disso, utilizaremos um nível de significância de 0.05 para os testes de hipóteses, o que nos dá um quantil de 1.96 para a distribuição normal padrão.
A simulação dos dados será feita a partir da biblioteca BTSR: Bounded Time Series Regression.
nH0 <- 100
nH1 <- 200
phi_parametro <- 0.2
alpha_teste <- 0.05
quantil_teste <- qnorm(1 - alpha_teste / 2)
coeficientes <- function(phi) {
# βARMA: the model from Rocha and Cribari-Neto (2009, 2017)
# is obtained by setting `coefs$d = 0`
# and `d = FALSE` and `error.scale = 1` (predictive scale)
list(alpha = 0, nu = 20, phi = phi, d = 0)
}
barma.sim <- function(n, phi, seed, y.start = NULL){
BARFIMA.sim(
n = n,
coefs = coeficientes(phi),
y.start = y.start,
error.scale = 1,
complete = F,
seed = seed
)
}
barma.phi_estimado <- function(yt, alpha = 0, nu = 20, phi = 0.1){
BARFIMA.fit(
yt = yt,
start = list(alpha = alpha, nu = nu, phi = phi),
p = 1, # Odem do polinômio AR
d = FALSE,
error.scale = 1,
report = F
)$coefficients["phi"][[1]]
}
barma.residuos <- function(yt, phi_estimado){
BARFIMA.extract(
yt = yt,
coefs = coeficientes(phi_estimado),
llk = F,
info = F,
error.scale = 1
)$error
}# Função para forçar a execução de `BARFIMA.extract` até que ela funcione
so_quero_que_funcione <- function(func, debug = F, max = 5, ...) {
FINALMENTE <- F
contador <- 0
while (!FINALMENTE) {
contador <- contador + 1
if (debug) print(paste("Tentativa:", contador))
resultado <- tryCatch({
func(...)
}, error = function(err) {
if (contador >= max) {
stop("Deu ruim, número máximo de tentativas atingido")
}
if (any(grepl("\\.btsr\\.extract\\(", deparse(err$trace$call)))) {
# Ignora erro de extração de resíduos
return(NULL)
}
stop(err)
})
if (!is.null(resultado)) {
FINALMENTE <- T
}
}
return(resultado)
}# Função auxiliar para gerar os dados
# Exemplo 1:
# ```r
# teste1 <- gerador_monte_carlo(
# parametros = list(
# lambda = seq(0.1, 0.4, by = 0.1)
# ),
# numero_de_execucoes = 2
# )
#
# teste1
# ```
#
# Exemplo 2:
# ```r
# teste2 <- gerador_monte_carlo(
# parametros = expand.grid(
# desvio_detectavel = seq(0.6, 1.8, by = 0.2),
# intervalo_de_decisao = seq(4, 6, by = 1)
# )
# numero_de_execucoes = 2
# )
#
# teste2
# ```
#
# Total de execuções é:
# `execucoes = numero_de_execucoes * length(novos_phis) * nrow(parametros)`
#
gerador_monte_carlo <- function(parametros = NULL, numero_de_execucoes = 200){
# Verifica se `parametros` é um data frame, caso contrário converte para um data frame
if (!is.null(parametros) && !is.data.frame(parametros)) {
parametros <- as.data.frame(parametros)
}
novos_phis <- seq(0.2, 0.6, by = 0.1)
# Expande a grade de parâmetros para cobrir todas as combinações
result <- expand.grid(
k = 1:numero_de_execucoes,
h1_phi = novos_phis
)
if (!is.null(parametros)) {
result <- merge(result, parametros, all = TRUE)
}
result %>%
mutate(id = row_number()) %>%
rowwise() %>%
mutate(
# H0
## Gera amostra de controle
h0_amostras = list(barma.sim(
n = nH0,
phi = phi_parametro,
seed = id
)),
## Estima o valor de phi da amostra de controle
h0_phi = barma.phi_estimado(h0_amostras),
## Calcula os resíduos da amostra de controle
h0_residuos = list(barma.residuos(h0_amostras, h0_phi)),
# H1
## Gera a amostra subsequente
h1_controle = list(barma.sim(
n = nH1,
phi = h1_phi,
# última observação da amostra de controle
y.start = h0_amostras[nH0],
seed = id + 1337E4
)),
## Calcula os resíduos da amostra subsequente
h1_residuos = list(so_quero_que_funcione(
\() barma.residuos(h1_controle, h0_phi)
))
)
}Vamos verificar nossa implementação do gerador de Monte-Carlo.
teste_montecarlo <- cache_dados(
"cache-teste-montecarlo",
function() {
gerador_monte_carlo(
numero_de_execucoes = 200
) %>%
select(h0_phi, h1_phi) %>%
mutate(
# Calcula a diferença entre os valores de phi
# `phi_parametro`: valor de phi definido para a amostra de controle
# `h0_phi`: valor de phi estimado para a amostra de controle. Espera-se que seja igual a `phi_parametro`
diferenca = h0_phi - phi_parametro
)
}
)ggplotly(
teste_montecarlo %>%
ggplot(aes(x = factor(1), y = diferenca)) +
geom_boxplot() +
geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
annotate(
geom = "text",
x = 0.55,
y = 0 + 0.01,
label = "0.00",
color = "red"
) +
labs(
x = "Valores de Φ₀",
y = "Diferença entre Φ e Φ₀",
title = paste0("Diferença entre Φ e Φ₀, para ", nrow(teste_montecarlo), " execuções")
) +
theme.base
) %>%
plotly.baseO EWMA (Exponential Weighted Moving Average) é um método de controle de processos que utiliza uma média móvel ponderada exponencialmente para detectar mudanças no processo. Sendo definido por:
\[ z_i = \lambda y_i + (1 - \lambda) z_{i-1} \]
onde, \(z_i\) é a estatística de controle no instante \(i\), \(y_i\) é a observação no instante \(i\), \(\lambda\) é o fator de suavização e \(z_{i-1}\) é a estatística de controle no instante anterior.
O fator \(\lambda\) é uma constante definida no intervalo \((0, 1]\) e seu valor inicial (para \(i = 1\)) é definido como a média do processo, tal que \(z_0 = \mu_0\).
A estatística de controle, \(z_i\), é comparada com os limites de controle, \(\text{LCS}\) e \(\text{LCI}\), sendo é definido como:
\[ \text{LC} = \mu_0 \pm L \sigma \sqrt{\frac{\lambda}{2 - \lambda} \left[1 - \left(1 - \lambda\right)^{2i}\right]} \]
Aqui, utilizamos o pacote qcc para implementar o
EWMA.
ewma_qcc <- function(amostra_inicial, lambda, amostra_subsequente) {
ew <- ewma(
amostra_inicial,
lambda = lambda,
nsigmas = quantil_teste,
plot = F,
newdata = amostra_subsequente
)
registros <- (nH0 + 1):(nH0 + nH1)
ewma <- as.numeric(ew$y[registros])
LI <- ew$limits[, 1][registros]
LS <- ew$limits[, 2][registros]
fora_de_controle <- ewma < LI | ewma > LS
total_fora_de_controle <- sum(fora_de_controle, na.rm = TRUE)
fracao_fora_de_controle <- total_fora_de_controle / length(ewma)
list(
ewma = ewma,
LI = LI,
LS = LS,
fora_de_controle = fora_de_controle,
total_fora_de_controle = total_fora_de_controle,
fracao_fora_de_controle = fracao_fora_de_controle
)
}ewma_monte_carlo <- cache_dados(
"cache-simulacao-ewma",
function() {
gerador_monte_carlo(
parametros = list(
lambda = seq(0.1, 0.4, by = 0.1)
)
) %>%
mutate(
qcc = list(ewma_qcc(h0_residuos, lambda, h1_residuos)),
ewma = list(qcc$ewma),
LI = list(qcc$LI),
LS = list(qcc$LS),
fora_de_controle = list(qcc$fora_de_controle),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
) %>%
select(-qcc)
}
)λ = 0.2Vamos analisar o EWMA para \(\lambda = 0.2\). Aqui, vamos considerar apenas a primeira execução.
ggplotly(
ewma_monte_carlo %>%
filter(k == 1 & lambda == 0.2) %>% # Apenas a primeira execução
select(h1_phi, ewma, LI, LS, fora_de_controle, lambda) %>%
rename(`Φ₁` = h1_phi) %>%
mutate(
observacao = list(1:nH1)
) %>%
unnest(
cols = c(`Φ₁`, observacao, ewma, LI, LS, fora_de_controle)
) %>%
ggplot() +
geom_line(aes(x = observacao, y = LI, color = "Limite Inferior"), linetype = "dashed") +
geom_line(aes(x = observacao, y = LS, color = "Limite Superior"), linetype = "dashed") +
geom_line(aes(x = observacao, y = ewma, color = "EWMA")) +
geom_point(data = . %>% filter(fora_de_controle),
aes(x = observacao, y = ewma, color = "Fora de controle"), size = 1.5, shape = 4) +
labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
scale_color_manual(
values = c(
"Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
"Limite Superior" = adjustcolor("red", alpha.f = 0.5),
"EWMA" = adjustcolor("black", alpha.f = 0.8),
"Fora de controle" = adjustcolor("red", alpha.f = 0.4)
)
) +
labs(
title = "EWMA para resíduos βAR(1), com λ = 0.2",
) +
theme.base +
theme(
legend.position = "bottom",
legend.title = element_blank()
) +
facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
plotly.baseResumo da fração de pontos fora de controle para diferentes valores de \(\lambda\) e \(\Phi_1\).
ewma_monte_carlo_resumo <- ewma_monte_carlo %>%
group_by(lambda, h1_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
datatable(
ewma_monte_carlo_resumo %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "Pontos fora de controle por Φ₁ e λ",
colnames = c("Valores de λ", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)Análise gráfica das médias de FPFCs (Fração de Pontos Fora de Controle).
ggplotly(
ewma_monte_carlo_resumo %>%
ggplot(aes(x = h1_phi, y = mean, color = factor(lambda))) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
color = "Valores de λ",
title = "FPFCs por Φ e λ"
) +
theme.base
) %>%
plotly.baseE, a distribuiçõa dos PFCs.
ggplotly(
ewma_monte_carlo %>%
ggplot(aes(x = factor(h1_phi), y = fracao_fora_de_controle, fill = factor(lambda))) +
geom_boxplot() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
fill = "Valores de λ",
title = "FPFCs por Φ e λ"
) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')O CUSUM (Cumulative Sum) é um método de controle de processos que utiliza a soma acumulada dos resíduos para detectar mudanças no processo.
Neste método são definidas duas constantes, \(H\) e \(K\), que representam o tamanho da mudança que se deseja detectar e a sensibilidade do método, respectivamente. As estatísticas de controle, sendo duas, são definidas como:
\[ \begin{matrix} \text{C}^+_i & = & \max\left[0, x_i - (\mu_0 + K) + C^+_{i-1}\right] \\ \text{C}^-_i & = & \max\left[0, (\mu_0 - K) - x_i + C^-_{i-1}\right] \\ \end{matrix} \]
onde \(x_i\) é a observação no instante \(i\), \(\mu_0\) é a média do processo, \(C^+_0 = C^-_0 = 0\) e \(i = 1, 2, \ldots, n\).
A constante \(K\), chamada de valor de referência, geralmente, segundo Montgomery (2009), é definida como a metade entre o \(\mu_0\) alvo e o valor fora de controle de média \(\mu_1\) que queremos detectar.
Já a constante \(H\), chamada de limite de decisão, é definida como o valor que a estatística de controle deve atingir para que possamos detectar a mudança no processo. Para Montgomery (2009), um valor razoável para \(H\) é \(5\sigma\). Assim, se \(\text{C}^+_i \geq H\) ou \(\text{C}^-_i \leq -H\), então o processo está fora de controle.
cusum_qcc <- function(amostra_inicial_, desvio_detectavel = 1, intervalo_de_decisao = 5, amostra_subsequente = NULL) {
arguments <- list(
data = amostra_inicial_,
se.shift = desvio_detectavel,
decision.interval = intervalo_de_decisao,
plot = F
)
if (!is.null(amostra_subsequente)) {
arguments$newdata <- amostra_subsequente
}
cu <- do.call(cusum, arguments)
registros <- if (is.null(amostra_subsequente)) {
seq(1, nH0)
} else {
seq(nH0 + 1, nH0 + nH1)
}
pos <- cu[["pos"]][registros]
neg <- cu[["neg"]][registros]
fora_de_controle_pos <- pos > intervalo_de_decisao
fora_de_controle_neg <- neg < -intervalo_de_decisao
total_fora_de_controle <- sum(fora_de_controle_pos) + sum(fora_de_controle_neg)
fracao_fora_de_controle <- total_fora_de_controle / length(registros)
list(
pos = pos,
neg = neg,
fora_de_controle_pos = fora_de_controle_pos,
fora_de_controle_neg = fora_de_controle_neg,
total_fora_de_controle = total_fora_de_controle,
fracao_fora_de_controle = fracao_fora_de_controle
)
}KQueremos que, sob \(H_1\), nas 100
amostras iniciais a probabilidade de um ponto fora de controle seja de,
nominalmente, 5%. Portanto, vamos buscar um valor para K
onde isso acontece.
desvio_detectavel_otimo_optimize <- function(amostra_inicial) {
for (x in seq(0, 3, by = 0.1)) {
fora <- cusum_qcc(amostra_inicial, desvio_detectavel = x)$total_fora_de_controle
if (fora == 5) {
return(x)
}
}
return(NA)
}
cusum_k_otimo <- cache_dados(
"cache-cusum-k-otimo",
\() {
data.frame(
id = 1:1000
) %>%
rowwise() %>%
mutate(
h0_amostra = list(barma.sim(
n = nH0,
phi = phi_parametro,
seed = id
)),
h0_phi = list(barma.phi_estimado(h0_amostra)),
h0_residuos = list(barma.residuos(h0_amostra, h0_phi)),
minimo = desvio_detectavel_otimo_optimize(h0_residuos)
)
}
)ggplotly(
cusum_k_otimo %>%
filter(!is.na(minimo)) %>%
ggplot(
aes(x = 0, y = minimo, group = 0, fill = "red")
) +
geom_boxplot() +
labs(
x = element_blank(),
y = "Valor ótimo de K",
title = "Valor ótimo de K por simulação"
) +
theme.base +
theme.no_legend
) %>%
plotly.basevalor_otimo_k <- mean(cusum_k_otimo$minimo, na.rm = TRUE)
kableExtra::kbl(
cusum_k_otimo %>%
filter(!is.na(minimo)) %>%
group_by() %>%
summarise(
"Média" = mean(minimo),
"Desvio padrão" = sd(minimo),
"Amostras" = n(),
),
caption = "Valor ótimo para K",
align = 'c',
digits = 3
)| Média | Desvio padrão | Amostras |
|---|---|---|
| 0.507 | 0.193 | 260 |
cumsum_k_monte_carlo <- cache_dados(
"cache-simulacao-cusum-k",
function() {
gerador_monte_carlo(
parametros = list(
desvio_detectavel = seq(0.1, 1.0, by = 0.1)
)
) %>%
mutate(
qcc = list(cusum_qcc(h0_residuos, desvio_detectavel = desvio_detectavel, amostra_subsequente = h1_residuos)),
pos = list(qcc$pos),
neg = list(qcc$neg),
fora_de_controle_pos = list(qcc$fora_de_controle_pos),
fora_de_controle_neg = list(qcc$fora_de_controle_neg),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
) %>%
select(-qcc)
}
)K ótimoVamos analisar o CUSUM para o valor ótimo de K
encontrado.
cusum_para_k_otimo <- cumsum_k_monte_carlo %>%
filter(k == 1 & desvio_detectavel == round(valor_otimo_k, 1)) %>%
rename(`Φ₁` = h1_phi) %>%
mutate(
observacao = list(1:nH1)
)
ggplotly(
cusum_para_k_otimo %>%
unnest(
cols = c(`Φ₁`, observacao, pos, neg, h1_residuos, fora_de_controle_pos, fora_de_controle_neg)
) %>%
ggplot() +
geom_hline(aes(yintercept = 5, color = "Limite de decisão"), linetype = "dashed") +
geom_hline(aes(yintercept = -5, color = "Limite de decisão"), linetype = "dashed") +
geom_line(aes(x = observacao, y = pos, color = "N+")) +
geom_line(aes(x = observacao, y = neg, color = "N-")) +
geom_point(data = . %>% filter(fora_de_controle_pos),
aes(x = observacao, y = pos, color = "Fora de controle"), size = 1.5, shape = 4) +
geom_point(data = . %>% filter(fora_de_controle_neg),
aes(x = observacao, y = neg, color = "Fora de controle"), size = 1.5, shape = 4) +
labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
scale_color_manual(
values = c(
"N+" = adjustcolor("blue", alpha.f = 0.6),
"N-" = adjustcolor("darkgreen", alpha.f = 0.6),
"Fora de controle" = adjustcolor("red", alpha.f = 0.3),
"Fora de controle" = adjustcolor("red", alpha.f = 0.4),
"Limite de decisão" = adjustcolor("red", alpha.f = 0.5)
)
) +
labs(
title = "CUSUM para resíduos βAR(1)",
) +
theme.base +
theme(
legend.position = "bottom",
legend.title = element_blank()
) +
facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
plotly.baseEm média, diminuir o K aumenta a sensibilidade do CUSUM
para detectar mudanças no processo.
cumsum_k_monte_carlo_resumo <- cumsum_k_monte_carlo %>%
group_by(desvio_detectavel, h1_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
datatable(
cumsum_k_monte_carlo_resumo %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "Pontos fora de controle por Φ₁ e K, com H = 5",
colnames = c("Valores de K", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)Análise gráfica das médias de FPFCs.
ggplotly(
cumsum_k_monte_carlo_resumo %>%
ggplot(aes(x = h1_phi, y = mean, color = factor(desvio_detectavel))) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
color = "Valores de K",
title = "FPFCs por Φ₁ e K com H = 5"
) +
theme.base
) %>%
plotly.baseE, a distribuição dos PFCs.
ggplotly(
cumsum_k_monte_carlo %>%
ggplot(aes(x = factor(h1_phi), y = fracao_fora_de_controle, fill = factor(desvio_detectavel))) +
geom_boxplot() +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
fill = "Valores de K",
title = "FPFCs por Φ₁ e K com H = 5"
) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')HSemelhante ao valor ótimo para K, queremos que, sob
\(H_1\), nas 100 amostras iniciais a
probabilidade de um ponto fora de controle seja de, nominalmente, 5%.
Portanto, vamos buscar um valor para H onde isso
acontece.
intervalo_de_decisao_otimo_optimize <- function(amostra_inicial) {
for (x in seq(1, 6, by = 1)) {
fora <- cusum_qcc(amostra_inicial, intervalo_de_decisao = x)$total_fora_de_controle
if (fora == 5) {
return(x)
}
}
return(NA)
}
cusum_h_otimo <- cache_dados(
"cache-cusum-h-otimo",
\() {
data.frame(
id = 1:1000
) %>%
rowwise() %>%
mutate(
h0_amostra = list(barma.sim(
n = nH0,
phi = phi_parametro,
seed = id
)),
h0_phi = list(barma.phi_estimado(h0_amostra)),
h0_residuos = list(barma.residuos(h0_amostra, h0_phi)),
minimo = intervalo_de_decisao_otimo_optimize(h0_residuos)
)
}
)ggplotly(
cusum_h_otimo %>%
filter(!is.na(minimo)) %>%
ggplot(
aes(x = 0, y = minimo, group = 0, fill = "red")
) +
geom_boxplot() +
labs(
x = element_blank(),
y = "Valor ótimo de K",
title = "Valor ótimo de K por simulação"
) +
theme.base +
theme.no_legend
) %>%
plotly.basevalor_otimo_h <- mean(cusum_h_otimo$minimo, na.rm = TRUE)
kableExtra::kbl(
cusum_h_otimo %>%
filter(!is.na(minimo)) %>%
group_by() %>%
summarise(
"Média" = mean(minimo),
"Desvio padrão" = sd(minimo),
"Amostras" = n(),
),
caption = "Valor ótimo para H",
align = 'c',
digits = 3
)| Média | Desvio padrão | Amostras |
|---|---|---|
| 3.16 | 0.712 | 125 |
cumsum_h_monte_carlo <- cache_dados(
"cache-simulacao-cusum-h",
function() {
gerador_monte_carlo(
parametros = list(
intervalo_de_decisao = seq(1, 6, by = 1)
)
) %>%
mutate(
qcc = list(cusum_qcc(h0_residuos, intervalo_de_decisao = intervalo_de_decisao, amostra_subsequente = h1_residuos)),
pos = list(qcc$pos),
neg = list(qcc$neg),
fora_de_controle_pos = list(qcc$fora_de_controle_pos),
fora_de_controle_neg = list(qcc$fora_de_controle_neg),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
) %>%
select(-qcc)
}
)K ótimoVamos analisar o CUSUM para o valor ótimo de K
encontrado.
cusum_para_h_otimo <- cumsum_h_monte_carlo %>%
filter(k == 1 & intervalo_de_decisao == round(valor_otimo_h)) %>%
rename(`Φ₁` = h1_phi) %>%
mutate(
observacao = list(1:nH1)
)
ggplotly(
cusum_para_h_otimo %>%
unnest(
cols = c(`Φ₁`, observacao, pos, neg, h1_residuos, fora_de_controle_pos, fora_de_controle_neg)
) %>%
ggplot() +
geom_hline(aes(yintercept = 5, color = "Limite de decisão"), linetype = "dashed") +
geom_hline(aes(yintercept = -5, color = "Limite de decisão"), linetype = "dashed") +
geom_line(aes(x = observacao, y = pos, color = "N+")) +
geom_line(aes(x = observacao, y = neg, color = "N-")) +
geom_point(data = . %>% filter(fora_de_controle_pos),
aes(x = observacao, y = pos, color = "Fora de controle"), size = 1.5, shape = 4) +
geom_point(data = . %>% filter(fora_de_controle_neg),
aes(x = observacao, y = neg, color = "Fora de controle"), size = 1.5, shape = 4) +
labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
scale_color_manual(
values = c(
"N+" = adjustcolor("blue", alpha.f = 0.6),
"N-" = adjustcolor("darkgreen", alpha.f = 0.6),
"Fora de controle" = adjustcolor("red", alpha.f = 0.3),
"Fora de controle" = adjustcolor("red", alpha.f = 0.4),
"Limite de decisão" = adjustcolor("red", alpha.f = 0.5)
)
) +
labs(
title = "CUSUM para resíduos βAR(1)",
) +
theme.base +
theme(
legend.position = "bottom",
legend.title = element_blank()
) +
facet_wrap(vars(`Φ₁`), ncol = 2, labeller = "label_both")
) %>%
plotly.baseEm média, diminuir o H aumenta a sensibilidade do CUSUM
para detectar mudanças no processo.
cumsum_h_monte_carlo_resumo <- cumsum_h_monte_carlo %>%
group_by(intervalo_de_decisao, h1_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
datatable(
cumsum_h_monte_carlo_resumo %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "Pontos fora de controle por Φ₁ e H, com K = 1",
colnames = c("Valores de H", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)Análise gráfica das médias de FPFCs.
ggplotly(
cumsum_h_monte_carlo_resumo %>%
ggplot(aes(x = h1_phi, y = mean, color = factor(intervalo_de_decisao))) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
color = "Valores de H",
title = "FPFCs por Φ₁ e H, com K = 1",
) +
theme.base
) %>%
plotly.baseE, a distribuição dos PFCs.
ggplotly(
cumsum_h_monte_carlo %>%
ggplot(aes(x = factor(h1_phi), y = fracao_fora_de_controle, fill = factor(intervalo_de_decisao))) +
geom_boxplot() +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
fill = "Valores de K",
title = "FPFCs por Φ₁ e K com H = 5"
) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')Vamos comparar os algoritmos EWMA e CUSUM(K) e CUSUM(H) para diferentes valores de \(\Phi_1\).
Aqui chamamos de CUSUM(K) as estatísticas onde o \(H\) for mantido constante e o \(K\) é variado, e o inverso para CUSUM(H).
estatisticas_resumo_combinadas <- bind_rows(
cumsum_k_monte_carlo_resumo %>%
filter(desvio_detectavel == 0.9) %>%
rename(parametro = desvio_detectavel) %>% mutate(algoritmo = "CUSUM(K)"),
cumsum_h_monte_carlo_resumo %>%
filter(intervalo_de_decisao == 4) %>%
rename(parametro = intervalo_de_decisao) %>% mutate(algoritmo = "CUSUM(H)"),
ewma_monte_carlo_resumo %>%
filter(lambda == 0.4) %>%
rename(parametro = lambda) %>% mutate(algoritmo = "EWMA")
) %>%
select(algoritmo, parametro, h1_phi, mean) %>%
mutate(
parametro = as.factor(parametro)
)ggplotly(
estatisticas_resumo_combinadas %>%
ggplot() +
geom_line(
aes(x = h1_phi, y = mean, color = parametro, linetype = algoritmo)
) +
geom_point(
aes(x = h1_phi, y = mean, color = parametro, shape = algoritmo)
) +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
linetype = "Algoritmo",
title = "Fração de pontos fora de controle",
color = "Valores dos parâmetros",
shape = "Algoritmo"
) +
theme.base
) %>%
plotly.baseestatisticas_combinadas <- bind_rows(
cumsum_k_monte_carlo %>%
filter(desvio_detectavel == 0.9) %>%
rename(parametro = desvio_detectavel) %>% mutate(algoritmo = "CUSUM(K)"),
cumsum_h_monte_carlo %>%
filter(intervalo_de_decisao == 4) %>%
rename(parametro = intervalo_de_decisao) %>% mutate(algoritmo = "CUSUM(H)"),
ewma_monte_carlo %>%
filter(lambda == 0.4) %>%
rename(parametro = lambda) %>% mutate(algoritmo = "EWMA")
) %>%
select(algoritmo, parametro, h1_phi, fracao_fora_de_controle) %>%
mutate(
parametro = as.factor(parametro),
h1_phi = as.factor(h1_phi)
)ggplotly(
estatisticas_combinadas %>%
ggplot() +
geom_boxplot(
aes(x = h1_phi, y = fracao_fora_de_controle, fill = parametro, shape = algoritmo),
) +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
fill = "Valores dos parâmetros",
linetype = "Algoritmo",
title = "Fração de pontos fora de controle",
color = "Valores dos parâmetros"
) +
facet_wrap(~algoritmo) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')A partir da análise da média dos pontos fora de controle, observa-se que o CUSUM(K), em média, possui poder semelhante ao CUSUM(H), enquanto que o EWMA é inferior aos dois.
Já a análise gráfica dos boxplots mostra que o CUSUM(K) possui uma variabilidade, em geral, maior que o CUSUM(H), enquanto que o EWMA possui, em geral, a menor variabilidade entre os algoritmos.
Vamos verificar se existe alguma combinação ótima entre os valores de
K e H para detectar pontos fora de
controle.
comb_cusum <- cache_dados(
"cache-simulacao-cusum-k-h",
function() {
gerador_monte_carlo(
parametros = expand.grid(
desvio_detectavel = seq(0.6, 1.8, by = 0.2),
intervalo_de_decisao = seq(3, 6, by = 1)
)
) %>%
mutate(
qcc = list(cusum_qcc(h0_residuos,
desvio_detectavel = desvio_detectavel,
intervalo_de_decisao = intervalo_de_decisao,
amostra_subsequente = h1_residuos)),
pos = list(qcc$pos),
neg = list(qcc$neg),
fora_de_controle_pos = list(qcc$fora_de_controle_pos),
fora_de_controle_neg = list(qcc$fora_de_controle_neg),
total_fora_de_controle = qcc$total_fora_de_controle,
fracao_fora_de_controle = qcc$fracao_fora_de_controle
) %>%
select(-qcc)
}
)comb_cusum_resumo <- comb_cusum %>%
group_by(desvio_detectavel, intervalo_de_decisao, h1_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
) %>%
mutate(
desvio_detectavel = as.factor(desvio_detectavel),
intervalo_de_decisao = as.factor(intervalo_de_decisao),
)
datatable(
comb_cusum_resumo %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "Pontos fora de controle",
colnames = c("Valores de K", "Valores de H", "Valores de Φ₁", "Média", "Mínimo", "Máximo")
)ggplotly(
comb_cusum_resumo %>%
group_by(desvio_detectavel, intervalo_de_decisao) %>%
summarise(
h1_phi = list(h1_phi),
mean = list(mean),
.groups = "drop"
) %>%
filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.06)) %>%
unnest(cols = c(h1_phi, mean)) %>%
ggplot(aes(x = h1_phi, y = mean, color = desvio_detectavel, linetype = intervalo_de_decisao)) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
color = "Valores de K",
linetype = "Valores de H",
title = "FPFCs por Φ₁, K e H"
) +
theme.base
) %>%
plotly.baseDe uma forma geral, observa-se que aumentando K e
H diminui a sensibilidade do CUSUM para detectar pontos
fora de controle.
E, analisando juntamente com as outras execuções, temos
estatisticas_combinadas_resumo_kh <- bind_rows(
estatisticas_resumo_combinadas,
comb_cusum_resumo %>%
filter(desvio_detectavel == 0.8 & intervalo_de_decisao == 5) %>%
mutate(parametro = as.factor(paste0(desvio_detectavel, ";", intervalo_de_decisao))) %>%
mutate(algoritmo = "CUSUM(K/H)") %>%
select(algoritmo, parametro, h1_phi, mean)
)ggplotly(
estatisticas_combinadas_resumo_kh %>%
ggplot() +
geom_line(
aes(x = h1_phi, y = mean, color = parametro, linetype = algoritmo)
) +
geom_point(
aes(x = h1_phi, y = mean, color = parametro, shape = algoritmo)
) +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
linetype = "Algoritmo",
title = "Fração de pontos fora de controle",
color = "Valores dos parâmetros",
shape = "Algoritmo"
) +
theme.base
) %>%
plotly.baseestatisticas_combinadas_kh <- bind_rows(
estatisticas_combinadas,
comb_cusum %>%
filter(desvio_detectavel == 0.8 & intervalo_de_decisao == 5) %>%
mutate(parametro = paste0(desvio_detectavel, ";", intervalo_de_decisao)) %>%
mutate(algoritmo = "CUSUM(K/H)") %>%
select(algoritmo, parametro, h1_phi, fracao_fora_de_controle) %>%
mutate(
parametro = as.factor(parametro),
h1_phi = as.factor(h1_phi)
)
)ggplotly(
estatisticas_combinadas_kh %>%
ggplot() +
geom_boxplot(
aes(x = h1_phi, y = fracao_fora_de_controle, fill = parametro, shape = algoritmo),
) +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
fill = "Valores dos parâmetros",
linetype = "Algoritmo",
title = "Fração de pontos fora de controle",
color = "Valores dos parâmetros"
) +
facet_wrap(~algoritmo) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')De uma forma geral, observa-se um leve aumento no poder quando
combinamos H = 5 e K = 0.8. Além diss, houve
um leve aumento na variabilidade.
Vamos combinar os dois algoritmos e ver o que acontece. A ideia é que, se um dos algoritmos detectar um ponto fora de controle, então o processo está fora de controle.
comb_monte_carlo <- cache_dados(
"cache-simulacao-ewma-cusum",
function() {
gerador_monte_carlo(
parametros = expand.grid(
desvio_detectavel = seq(1.6, 2.0, by = 0.2),
intervalo_de_decisao = seq(4, 6, by = 1),
lambda = seq(0.6, 1.0, by = 0.2)
)
) %>%
mutate(
"(K;H;λ)" = factor(paste0(desvio_detectavel, ";", intervalo_de_decisao, ";", lambda)),
ewma = list(ewma_qcc(h0_residuos, lambda, h1_residuos)),
cumsum = list(cusum_qcc(h0_residuos,
desvio_detectavel = desvio_detectavel,
intervalo_de_decisao = intervalo_de_decisao,
amostra_subsequente = h1_residuos)),
fora_de_controle = list(
ewma$fora_de_controle | (cumsum$fora_de_controle_pos | cumsum$fora_de_controle_neg)
),
fracao_fora_de_controle = sum(fora_de_controle) / nH1
) %>%
select(-ewma, -cumsum)
}
)comb_monte_carlo_resumo <- comb_monte_carlo %>%
group_by(`(K;H;λ)`, h1_phi, desvio_detectavel, lambda) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
datatable(
comb_monte_carlo_resumo %>%
select(-`(K;H;λ)`) %>%
mutate(across(where(is.numeric), \(x) round(x, digits = 4))),
caption = "Pontos fora de controle por Φ₁ e (K;H;λ)",
colnames = c("Valores de Φ₁", "K", "H", "λ", "Média", "Mínimo", "Máximo")
)ggplotly(
comb_monte_carlo_resumo %>%
# get `(K;H;λ)` where h1_phi == 0.2 and mean >= 0.05 & mean < 0.06
group_by(`(K;H;λ)`) %>%
summarise(
h1_phi = list(h1_phi),
mean = list(mean),
.groups = "drop"
) %>%
filter(map_lgl(mean, ~ .x[1] >= 0.05 & .x[1] < 0.06)) %>%
unnest(cols = c(h1_phi, mean)) %>%
ggplot(aes(x = h1_phi, y = mean, color = `(K;H;λ)`)) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
color = "Constantes (K;H;λ)",
title = "FPFC por Φ₁ e (K;H;λ)"
) +
theme.base
) %>%
plotly.baseestatisticas_combinadas_resumo_ew_sum <- bind_rows(
estatisticas_combinadas_resumo_kh,
comb_monte_carlo_resumo %>%
filter(`(K;H;λ)` == "1.8;6;0.8") %>%
mutate(parametro = as.factor(`(K;H;λ)`)) %>%
mutate(algoritmo = "EW-SUM") %>%
select(algoritmo, parametro, h1_phi, mean)
)ggplotly(
estatisticas_combinadas_resumo_ew_sum %>%
ggplot() +
geom_line(
aes(x = h1_phi, y = mean, color = parametro, linetype = algoritmo)
) +
geom_point(
aes(x = h1_phi, y = mean, color = parametro, shape = algoritmo)
) +
labs(
x = "Valores de Φ₁",
y = "Fração de pontos fora de controle",
linetype = "Algoritmo",
title = "Fração de pontos fora de controle",
color = "Valores dos parâmetros",
shape = "Algoritmo"
) +
theme.base
) %>%
plotly.baseDe uma forma geral, houve uma piora na detecção de pontos fora de controle ao combinar os dois algoritmos. Isso pode ser explicado pelo fato de que, a sensibilidade do EWMA é maior que a do CUSUM, o que pode ter levado a um aumento de falsos positivos quando o processo permanece sob controle.